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ABSTRACT 


In this thesis we present an algorithm for the precise determination of Fermi-Dirac 
(FD) integral functions, f,(n), for arbitrary values of the parameter t and the argument n. 
The FD integrals are a class of functions that are used extensively in the modeling of 
semiconductor devices, e.g., when the charge carriers are in a strongly quantum, 
degenerate regime, such as in heavily doped semiconductors. The determination of FD 
integrals has a long history. Our approach to evaluating these functions is two-fold. First, 
we develop exact power series expansions of the integral. These series, however, 
converge too slowly to be a practical means of evaluating the integral. The second aspect 
of our approach 1s to apply numerical series acceleration methods to improve the rate of 
convergence of these power series expansions. Indeed, it would not be feasible to use the 
series expansions without implementing an acceleration method. The result is a computer 


program that provides efficient, accurate values of the FD integral. 
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I. INTRODUCTION 
In modeling semiconductor devices, the number of free electrons per unit volume, 
#7, in the conduction band depends on the density of states D(E) available for occupation, 
and on the Fermi Dirac distribution functionfpp(E), which yields the probability that a 


state of energy E is occupied [Ref. 1: pg. 90] 


E 


n= [dED(E) fpp(E) - (1) 


Here, E, is the conduction band energy lower level, and E; 1s the top of the conduction 


band. The three dimensional form of the density of states, D(E) is [Ref. 1: pg. 91] 


D(E) = = (22 © (ETE, (2) 








where m* is the effective mass of the charge carrier. The Fermi Dirac distribution function 
1S 


] 





ee 3 
J) = Ty exp(BE =H) ©) 
where f is equal to (kT)", and p1 is the Fermi energy. When E is equal to p, frp(E) is 
equal to 1/2. Combining Eqs. (1-3) 7 becomes 
) & ) 3/2 E,-E, dy.fy (4) 
n= i => == eee 
2) f+ exp(By +E. — 1) 


where y 1s equal to E-E,. 


Since the upper limit of the energy integral will represent few electrons, setting 


(E;-Ec) equal to infinity, 2, will introduce an exponentially small error into the value of the 
integral. This approximation 1s allowed because f(E) vanishes quickly for energy E >>. 


The equation for the number of carriers per unit volume is then 








| ——— dy/y 





R= so 5 
jew ijt , Repeat Eee) ©) 
We then have 
1 (2m*\ 

Mi a Bh F,.(B(u a E.)). (6) 

where Fy, denotes the following integral 

Zz 
© x 

Fi2(z)= I, ax . (7) 


]1+exp(x — z) 


This integral 1s but one example of an entire family of integrals labeled by the parameter t 


in the following manner 


t 


x 


ia (z) = i meme mek (8) 


These integrals are known as Fermi-Dirac integrals. The purpose of this thesis is to 
examine methods of computing Fermi-Dirac integrals. In the following sections we will 
use the conventional notation and take the definition of the F-D integral to be the 


following [Ref. 3: pg. 1067] 


| vee 


ea ae ° 
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where I is the standard gamma function [Ref. 2 . Ch. 6}. Returning to the case of the 


three dimensional electron gas, setting t =1/2; and using the fact '(3/2) becomes (11/2)”, 


we have 





3/2 
; | fir(B(E. - W)). (10) 


The equation for n becomes a constant times fy,(B(E.-)), as long as the temperature is 
invariant. 

Thus we see that the carrier density in a semiconductor is generally related to the 
FD integral function. Except 1n special cases, the FD integral cannot be evaluated 
analytically in closed form. In many treatments of the carrier density, the following 
approximation 1s adopted. When B(p-E.)<< -1, 1.e. when the Fermi level lies at least 
several multiples of kT below the conduction band, we can approximate f(y) ~ e’ for y 
<< 0. We will term this the Boltzmann limit, since in this case the Fermi Dirac distribution 


goes over to the Boltzmann distribution function. We then recover the often used 


Pen a 


3/2 
(Bem) where N, is equal to = = . However, for accurate 
2 





approximation n =N,e ” 


modeling of the carrier density in semiconductor, or when the Fermi level lies close to the 
conduction band edge, we can no longer invoke this approximation. Instead we must 
evaluate the Fermi-Dirac integral. 

When the Fermi level, 1, is close to the conduction band, E., the electron gas 1s 
considered degenerate. If x, which corresponds to Bu, is much greater than one, then the 
evaluation of the Fermi-Dirac integral must be evaluated in the extreme quantum or strong 


degeneracy regime. Accurate evaluations for the FD integral in the extreme quantum 


region can be done with well known techniques, such as the Lhrenburg, Joyce-Dixon, 
Chang-Izabelle, or the Nilsson approximations [Ref. 1: pp. 111-114]. However, in the 
intermediate semi-classical regime, neither the Boltzmann or the strong degenerate 
approximations are applicable. 

The semi-classical regime is of particular importance to the field of semiconductor 
physics, where the controlled doping of the intrinsic material can range from the classic to 
the strong degenerate regime. Attempts to evaluate the FD in the past have included tables 
of values by McDougall and Stoner (1938) [Ref. 7] for the range -4< x < 20, in increments 
of Ax = 0.1. Entries were accurate to seven digits for t equal to 1/2 and 3/2, and five digit 
accuracy for t= -1/2. Rhodes (1950) [Ref. 8] provided formulas and tables of F,(x), for t 
=], 2,3, and 4 in the interval |x| < 4. 0. A parallel effort has been made to develop 
approximants which can replicate the values of the Fermi-Dirac integrals, f,(x) to 
exemplary accuracy for different ranges of x. Whereas tabular data is used in conjunction 
with interpolation methods; by incorporating an approximate formula within an algorithm, 
the values of F, can be calculated directly by the user. 

In terms of the analytical treatment of F,(x) the code developed by Dingle (1957) 
[Ref 9] is quite valuable. His work has provided the portion of an FD integral that was 
discarded in Sommerfield’s (1928) [Ref. 10] treatment of the extreme quantum regime. 
Especially effective approximants were given by Code and Thatcher (1967) [Ref. 11] for t 
= -1/2, 1/2, 3/2, and by Van Halen and Pulfrey (1985) [Ref. 12] for t = -1/2, 1/2, 312, 5/2, 
3, and 7/2. [Ref. 3: pg. 1071]. Blakemore has provided reference and tabulated material; 
[Ref. 4: Appendix B] and summarized much of the data on Fermi Dirac integrals through 


1981 [Ref 3] 


In view of the continuing need to generate accurate values for Fermi Dirac 
integrals for arbitrary values of x, we were motivated to reexamine this subject and to 
develop an efficient algorithm that delivers precise values. Our program can be employed 
to generate directly the values of desired FD integrals, or to generate tables of highly 
accurate values which can then be used in conjunction with standard interpolation 
methods. 

Our approach, rather than following the route of numerical integration adopted by 
McDougall and Stoner, consists of utilizing exact series expansions of f,(x) .The apparent 
difficulty that a given expansion converges very slowly is overcome by the powerful series 
acceleration method of Levin [Ref. 5] , that will be covered in Chapter III. To obtain very 
high accuracy in the value of a slowly converging series, it turns out to be sufficient to 
input into the Levin method only the first fifteen terms of a given series expansion. 

In Chapter II we summarize several known formulas for FD integrals which apply 
for the integer values of t. This is followed by the derivation of exact series expansion of 
an FD integral, which can be used in conjunction with the Levin series acceleration 
method for general values of t, as discussed in Chapter III. Without the aid of a series 
acceleration method, the various formulas we give would, from an practical view, be 
ineffective throughout the interval -1 <x < 10. This underscores the great usefulness of 
the Levin method for the evaluation of the general Fermi Dirac integral. Finally, Chapter 
IV summarizes our results. In the appendices we present our algorithm, Blakemore's 
tabulated data for the half integer values of t, and a comparison of our values to 


Blakemore's in the region of -4 < x < 4, for different values of half integer t. 


The necessity to understand the importance of Fermi-Dirac integrals 1s imperative 
in the modeling of semiconductor devices. For example in the determination of n, the 
number of charge carriers, t = 1/2 is used to model a three dimensional semiconductor. In 
the case of a two dimensional electron gas, found, for example, in the channel of a 
MOSFET device, the effective density of states must be changed to Nc = m*/(mh 7B) and 
t will be equal to 0. For a one dimensional electron gas, found in modern "quantum wire" 
devices, N. is equal to [m*/(1h 7B) ]"* and t is equal to - 1/2 [Ref. 1: Ch. 3]. We note 
that these values of t for a D-dimensional electron gas, t = (D-2)/2, arise from the 
assumption of a parabolic band structure function E(t) = h74°/2m*, appropriate to free 
particles. When one includes a more realistic, nonparabolic band structure, such as 1s 
required for PbTe, PbSe, and PbS; then n turns out to be given by a linear combination of 
fsn(m) and fsr(n) [Ref. 6: pg. 47]. In this thesis we include the values of f,(x) for t = 7/2, 
in order to have full comparison with Blakemore's tabulated data. In this age of 
miniaturization of semiconductor devices the applicability of the Boltzmann approximation 
will soon be exhausted. The ability to accurately model the number of free electrons will 
require the utilization of numerical algorithms. Also, 1n the future, 1f the need arises to find 
the Fermi-Dirac integrals of an arbitrary value of t, the program we present will be 


capable. 
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Hl. ANALYTICAL METHODS 

In the first part of this chapter we list several known, convenient formulas, for 
f(x) which apply for integer values of t. This 1s followed by exact series expansions for 
FD integrals for general t in the separate regime of positive and negative values of x. The 
expansion we give for x > 0 1s new. When used with the Levin method, these expansions 
can provide values which are correct to not less than six significant figures for any value of 
the argument. Throughout this thesis we restrict our attention to real values of x. 

For t = 0 one can evaluate the FD integral in closed form, fo(x) = In(Il+e’). Starting 
with this result and repeatedly invoking the relationship 0 f,/@x = f.1, one can derive 
explicit formulas for f,(x) for t equal to any negative integer. Formulas for t equal to a 
positive integer can be obtain from fo by repeated integrations, but the evaluation of these 
integrals in closed form can not be achieved. 


For negative values of x, for arbitrary t, there exists the convergent expansion 





00 Nad e 
f(x) = e™ , (x<0). (11) 
n=] n 


This result is arrived at by simply substituting the geometric series expansion into Eq. (9): 


a5 
e€ 





: ee ; (le) 
n=! 


lees 


which converges for all s > 0 as long as x < 0, and integrating term by term. For the 


special case x = 0 we have 
f.(0) = n(t+1), (13) 


where 1(s) is defined for any positive integer s by [Ref. 2: Ch. 23] 


ns) = (CI fmt (14) 


Because of the frequent occurrence of the quantity n(s) in the treatment of the general FD 
integral, we summarize its major properties. It 1s closely related to the Rieman zeta 


function, C(s) defined by 


Ae) lar f (15) 


m=| 


Specifically one has [Ref. 2: Ch. 23] 





n(s) -(1- = Jets). (16) 


Highly accurate values for the n(s) are given in Table 23.3 of reference 2 and can easily be 
reproduced by using the Levin method. The latter route was used within the ANSI C code 
used in our algorithm. 

It should be noted that the series on the right hand side of Eq. (14) converges only 
if the Re s > 0, as should be expected since the expansion given by Eq. (11) for f,(x) 1s 
similarly restricted. Even though the expansion in Eq. (11) converges for x < 0, for 
numerical purposes it 1s nominally useful only if x < -1, because the rate of convergence of 
the series decreases sharply as x approaches 0 from below. In practice this equation 1s 
usetul within the interval -1 < x < 0 only if the method of Levin is used . 

For x > 0 it is useful to rewrite f as the sum of two integrals, the first extending 
from 0 to x and the second extending from x to ©. In the first of these intervals we change 
of variable s = x(1-y) and use the identity I1/(e*+1 ) = 1-1/(e”+ 1). In the second integral 


we Change variable, s = x(I+y) This gives 


x7! 


Tg 





ee {1+(c +1)[4-(x)-B,Co)}, (17) 


where the functions A,(X) and B,(x) are defined by the relations 


Coes y)’ 


A:(x) = [, o~ ae 


(18) 


Game 
wa | 


B,(x) = Los (19) 


Note that the integral in Eq. (18) converges as long as x > 0, whereas the integral in 
equation (19) converges for positive and negative values of x. 

Before presenting our method for evaluating these functions we summarize their 
analytical properties and derive several useful expansion formulas. If t is a positive integer 
or zero, the function A,(x) has an isolated pole at x = 0 of degree t+]. This can be seen 
directly by making the change of variable u = xy. For all other non-integer values of 1, 
A,(x) has an isolated branch point at x = 0. The leading behavior for small positive x is 
obtained by noting that in this regime the convergence of the integral is achieved by the 
growth of the exponential term for values of y satisfying the inequality y >> 1/x, for 
which the numerator in equation (18) can be approximated by y ‘. Thus, for very small 
positive values of x we have the result A,(x) ~ Cx” | where C =T'(t + 1)n(t + 1). The 
full expansion for small x is fairly involved and we shall not give it here. 

A far more useful expansion is provided by the following. We expand e’/(1 + e”’) 
as a geometric series in powers of e” , which converges for all positive values of x, utilize 


Eqe (sales) iomeference 2: 


in dy(1+ y)'e® =U(1, 2+7, x), (20) 


where U denotes the standard irregular solution of the confluent hypergeometric equation. 
The major properties of this function are listed in Chapter 13 of reference 2. The final 


result 1S 


4,0) = J (-DEUC, 2 + 7, x). (21) 


k=1 
We have found that the simplest way to evaluate the U functions appearing in this 


equation accurately is to use the continued fraction representation 


Ud. 2 pyatt gO 7 aes ee (22) 
IZ 2 2 2 2 il ate 


This representation converges for all z. In fact, the rate of convergence 1s greater the 
larger the value of |z| If expanded in powers of I/z, the result is the standard asymptotic 
expansion, 1.e. Eq. (13 .5 .2) of reference 2 
UGlie2 + pe Deg OT Izy er” : (250 
Here (Z), 1s the Pochammer symbol defined by Zp =1, Z,= (z + 1)...(z2 +n- 1) forn=1, 2, 
_ Note that, whereas the asymptotic expansion diverges for all z, the continued fraction 
expansion, equation (22), converges for all z. 

Returning to the series in equation (21), for large values of the argument kx, the 
function U(1, 2 + t, kx) decreases to zero as I/(kx). As such, the series in question is only 
conditionally convergent. The direct method of summing the series as it stands is 
impractical, because if one were to sum the first N terms the resulting error in the estimate 
of the infinite series would only be of order of the last term retained, 1.e. O(//N). 


However, although this expansion is slowly convergent, it can be evaluated to very high 


accuracy using the Levin method, and then it 1s necessary to evaluate only the fifteen or so 
terms of the series in equation (21). 

We now discuss some of the properties of the function B,(x). As remarked earlier, 
the integral in Eq. (19) converges for both positive and negative real values of x. In fact, 
with the aid of the identity I/(e + 1) = 1 - I/(e + 1) one finds that 

Breas B-x) = W(t + 1). (24) 
The singularities of B,(x) closest to the origin in the complex x plane are the pair of branch 
points at t/z, as can be seen from the fact that for these choices of x the integrand in Eq. 
(19) has simple poles for y = 1. Hence B,(x) admits a convergent Maclaurin expansion 
within the circle of radius m centered about the origin. The form of the Maclaurin 
expansion 1s readily established. We substitute into Eq. (19) the expansion 


eee 1B x 25) 
go) ED nV (t+1+2n) © 





n=} 
which converges for real values of y in the interval {0, 1 ] as long as |x| <7. The 
quantities B2, are the Bernoulli numbers, referred to above. One thus arrives at the 
Maclaurin expansion 


ee io. ii oe 
nV (zt +1+2n) 





B(xy=4 TED” (26) 
n=\ 


(c +1) 2 
Although one can utilize this expansion along with the Levin method to obtain accurate 
values for values of x somewhat beyond x = 7, a different type of expansion is necessary 


tor larger x. In the following we provide a far more effective expansion. The derivation 


starts from the definition, Eq. (19), and we utilize Eq. (13.2.1) of reference 2, 


ell 


i; ary) ea LM, es » MC. pe a a (27) 


where M denotes the usual confluent hypergeometric function, 





oo a) 
M(a,c,z) = (2), Bans (28) 
2 k\(c), 
In this way we arrive at 
B(x)=[l/ (r+ py (—1)**" A412 + c.—hex). (29) 


k=) 
Upon utilizing the asymptotic property [Ref. 2: Eq. (13.1.5) ] 

M(a, c, z) =[I'(c)/I'(c-a)](-z)“[1+O(1/z)] (Re z < 0), (30) 
one notes that the typical term of the series, including the factor (-1) Kl decreases to zero 
for large values of k proportional to (-l)‘" /k. It thus follows that the series in Eq.(29) is 
conditionally convergent for all positive values of x. As in the case of the series, Eq.(21), 
for A,(x), the Levin method can be used with great effectiveness for this conditionally 
convergent series. 

To evaluate the confluent hypergeometric functions M(1, 2 + t, -kx) appearing in 
Eq.(26) accurately, two routes are open. One method consists of exploiting the Kummer 
idemuty (KRete2- Eq. (isan 27 

M(a, c, z) = e’M(c-a, c, -z), (31) 


so that evaluation of B,(x) would proceed by invoking the Levin method on the series 
Bx) ={ WctH1)] So (-*e M+ 7,24+7,kx), (32) 
ko] 


A second. more effective, method for evaluating M is to use the continued fraction 


representation 


M(1, ¢, a ca, (33) 


ies jet+2 |ct+3 [e+4 
which converges for all finite z. Once accurate values of the first fifteen terms of the series 
in Eq. (32) are computed, the application of the Levin method to these terms of the series 
provides estimates of B,(x) of very high accuracy. 

In the following chapter we digress to provide a brief description of the Levin 
method and the conditions to be met in order for the method to be effective. We will show 
that these conditions are met for all of the expansions encountered in this work. In 
Chapter IV we display some numerical results obtained using the present method and 


discuss the accuracy that can be achieved. 
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I. SERIES ACCELERATION METHOD 

In this chapter we shall present a brief description of the Levin series acceleration 
method, the key formulas, as well as a list of conditions to be met for the method to be 
effective Although, the method first appeared tn print over two decades ago, it is still 
virtually unknown in the physics community. This method is for many classes of problems 
perhaps the most effective technique currently available to accelerate the convergence of a 
series. For further details on the general method the reader can consult Levin's original 
paper [Ref. 5]. 

The central goal is, given a quantity 7 represented by an infinite series, to obtain 
highly accurate estimates of the value of that quantity utilizing only the first few terms of 
the series. To fix the notation, let the individual terms of the infinite series be denoted by t, 


so that 


We will require values of the partial sums, T, (n = 1,2,3,...), of the infinite series, which 


are defined by 
va (35) 
k=) 


In the following, we explicitly consider the vast class of infinite series with the property 
that the error in approximating T by the value of the partial sum T, , is of order the last 
individual term retained, t,. As discussed below, this class includes convergent as well as 
divergent series. A quantitative statement of this property is provided by the relation 


- | ail a oak a (36) 


Ie 


for all positive integers n, where g, possesses a unique, finite limit for n—90. Now the 
sequence {g,} will possess this property as long as the ratios of successive terms of the 
given series, namely tp/ta+1 , approach a unique finite limit, differing from unity, for n>. 


To prove this, note that from Eq. (36) we have 


n+) 


f 
Sat = eae eZ (37) 
We now suppose that the ratio in parenthesis in Eq (37) has a unique limit, €, such that 
e = lim (aa) (3 8) 


exists. It thus follows from Eq. (37) and (38) that 


| 


(39) 


Lo— lim g, = 


n—>O l- 


Uy 


We therefore arrive at the claimed result that, if & exists and 1s unique and is unequal to 
unity, the sequence {g,} possesses a unique limit for n — ©, given by Eq. (39). 

It is important to note that Eqs. (37) and (38) can apply even if the infinite series 
representation of T diverges. For example, for the terms of the infinite series in Eq (11), 
which converge only if x < 0 one has & = -e™ so that g, =I/(I+e”). In particular, the 
quantity g,. exists and ts finite for both positive and negative values of x. Similarly, the 
infinite series representation of T can be a divergent asymptotic series, arrived at, for 
example, by formal expansion of a portion of the integrand of a bona fide integral 
representation of T. 

Now, the Levin method consists of adopting the assumption that g, is well 
approximated by a polynomial in I/n. We define a series of approximations g,[N] to the 


truc gy, defined in Eq. (36), by the following, 


N-2 


CHUN | = ee, (I/n)" , (40) 


k=0 
where N denotes the Nth Levin approximant. Inspecting Eq. (40) one can anticipate that 
this basic assumption of the Levin method will be appropriate if t,/t,.; is a slowly varying 
function of 1/n. Adopting Eq. (36) we proceed to require that for a fixed value of N, T, = 
T[N] + g,[N] for n= 1, 2, ....,N. Here T[N] denotes the Levin estimates of T as obtained 
by employing Eq. (36) and the first N terms of the series of Eq. (34) . We now have N 
equations which determine the N unknown quantities T[N], {cx}, (k = 0, 1, 2, .., N-2). It 


turns out that T[N] is given by 


T[N] = PINVQIN] (41) 
Where 
BUNS Senter, /t, NU[KICN -k)!] (42) 
= 
and 
Q[N] = Set N URI DY (43) 


k= 

We reiterate, the Levin estimates for the value of the infinite series T, which utilize 
the first N terms of the series, is given by equations (41-43). This estimate can be expected 
to serve as a good approximation for T in those cases where the ratios of successful terms, 
tp/tn+1, 18 a Slowly varying function of I/n. In this regard for the series expansions of an FD 
integral which have appeared in Chapter II, ( Eqs. 11, 21, 26, 29, and 32) the ratios of 
successive terms are in fact slowly varying in I/n. 

In practice where deployment of the method is justified, the estimates T[N] for N 


equal to 13, 14, or 15 typically are very nearly equal, and the first set of common digits 


can be taken as providing a reliable estimate of the value of the infinite series T. It is 
frequently impractical to consider larger values of N because the individual terms of the 
series in Eqs. (42) and (43) become so large the severe round off problems arise in the 
course of summing the series. 

For example, the value n(1) is In(2) or 0.693 147805599453 [Ref. 2]. By using Eq. 


(14) and summing the first 15 terms the value of n(1) is quickly resolved. As can be seen 







Se 01 
- 1 1-000000 | 1-000000 | 70000000000000000— 
601616667 [0.166667 | 0.6931474019283138_ 
8 [0.634534 | -0.125000 | 0.6931471 800150043 _ 
9 [0.745635 [OTT | 0.6931471 806012292 


Table 1 





in the table | by the fifteenth term in the series there is agreement with the accepted 
value of In(2) out to the 16th decimal place. If one one directly summed the series 
represcntation of nQ), Eq. (14), it would take significantly more than fifteen terms to reach 
this level of accord. In fact we used the first billion terms in Eq. (14) and reached 
agreement only to the fourth decimal place. The exercise required two hours of CPU time 
on a SPARC 10 Sun workstation. Clearly using the series acceleration method of Levin is 


the proper course of action whcn summing these types of series. 
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IV. RESULTS 

The algorithm was built as a set of subroutines to be called by the parent equation. 
For example Eq. (17) requires a gamma function. In turn the Eqs. (18) and (19) that are 
called by Eq. (17) require hypergeometric and confluent hypergeometric functions. 
These subroutine were tested against known values in reference 2. Once the gamma, 
hypergeometric, confluent hypergeometric and eta functions were working perfectly they 
were combined to get the desired Eqs. (17), (18) and (19). 

Next we ran the algorithm from the interval -4 < x < 4 to replicate the tabulated 
data in Blakemore. The values of t that were of interest were the half integer values -1/2, 
1/2, 3/2, 5/2, and 7/2. In Figures 1- 5 we show the error of our code as compared to 
Blakemore's data. The error was calculated by the formula [(B-L)/B]*100 to get a 
percentage, where B is Blakemore's value and L is the output of our algorithm. It is 
obvious from inspection of the figures that the accord with Blakemore's data is to the last 
digit. We, of course, claim that our results are more accurate than published results of 
Blakemore; because of suspected round off error in his results. The flow of data, as a 
function of x, is shown in Figure 6, where the Blakemore values are plotted as straight 
lines and F(x), the output of our program is superimposed with circles. In order to fully 
represent the values of our data tabulated comparison to Blakemore's values are collected 
in Appendix C. It is quite clear that we had agreement with Blakemore to the last digit for 
f(x) for all the values of x . The strength of this program is the speed of the Levin 
method. In the fraction of a minute all the values of f(x) were generated between 


-4<x < 4. 
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Figure 1: Percent Error for t 
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Figure 2: Percent Error for t = 1/2 
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Percent Error for Tau = 3/2 


Percent Error (B-L)/B * 100 
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Figure 4: Percent Error for t = 5/2 
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Percent Error for Tau = 7/2 
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Figure 5: Percent Error for t = 7/2 

Prior to using the method of Levin we tried to duplicate the results by letting the 
infinite series run until they were in agreement in the sixteen decimal place. This proved 
to be taxing on both machine and programmer. Once the jobs were completed there was 
poor agreement with accept value as shown in Figure 7 for t = 1/2. This fluctuation about 
the accepted value of f(x) is due to round off error in summing the infinite series. Since 
the algorithm ran for days there was plenty of time and opportunity for errors to build up 
and truncate through the iteration. For other values of t the disparity between accepted 
and calculated values were just as poor for values of x > 0. Therefore 1n order for series 
approximations to be useful in the modeling of FD integral, the use of series accelerations 


methods such as Levin is paramount. 
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Blakemore vs. Levin for Half Integer Values of t 
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Figure 6: Comparison of Levin to Blakemore for half integer values of t 





Figure 7: Illustration of the error generated when the Levin Method ts not used 
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V. CONCLUSIONS 

In this thesis we have developed a new, highly efficient algorithm for evaluating 
Fermi-Dirac (FD) integrals, f(x), the definition of which is given in equation (9). 
Fermi-Dirac integrals are widely used in numerous aspects of semiconductor modeling, 
as we have discussed in Chapter I. Except in special cases, FD integrals cannot be given 
in closed form using elementary functions. There has been a long history of attempts to 
provide methods for approximately evaluating these integrals. The algorithm developed 
in this thesis is based on the following two-fold approach. 

First, we employ various series expansions, derived in Chapter II starting from 
the definition of the FD integral. For x < 0, for example, we employ the expansion given 
by equation (11). For x >0, however, we utilize equation (17), which in turn is evaluated 
using the series expansions given by equations (21) and (29). These expansions provide 
an alternate, exact representation of the Fermi-Dirac integral, each applicable in its 
respective domain of validity. While these expansions are mathematically exact, their 
utility is seemingly limited since their rate of convergence is typically so slow as to 
render them impractical as a means of evaluating FD integrals. 

The second part of our approach, therefore, is to use a novel series acceleration 
method as applied to the expansions derived in Chapter II. The general aspects of this 
method, due originally to Levin, Ref. [5], are discussed in Chapter III. An illustration of 
the power of the Levin method was given in Chapter III. There, we applied the series 
acceleration method to the evaluation of the eta function, defined by the power series 


given in equation (14). This expansion converges too slowly, to the point where 
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straightforward summation of (14) 1s futile, since, as the series is an alternating series, 
round-off errors ultimately limit computational accuracy. Summation of a billion terms 
in Eq. (14) produced only 4 digits of accuracy. Application of the Levin technique, 
however, readily produced 16 digit accuracy in only 16 iterations of the method. 

The results we obtain by applying the Levin method to the exact series expansions 
of the FD integrals are discussed in Chapter IV. We have compared our calculated 
results with those tabulated 1n the book by Blakemore, Ref. [4], which is a standard in the 
literature of semiconductor modeling. In Figure 6 we show Blakemore’s results in the 
solid curves; are our results indicated with circles. On the scale of Figure 6 it would 
appear there is perfect agreement. Blakemore, however, only cites values of FD integrals 
to four significant digits. In Figures 1-5, we illustrate percentage error of our results as 
compared with those of Blakemore. In all cases our results agree with those of 
Blakemore to the number of digits that Blakemore lists. 

We believe our work sets a new eae for the computation of FD integrals. 
We note that while the most useful FD integrals are those with the parameter t restricted 
to half-integer values (1.e., t = -1/2, 1/2, 3/2, etc.), our algorithm can evaluate FD 
integrals with arbitrary values of t, a feature not previously available to the best of our 
knowledge. From a practical point of view, one could envision two uses of the algorithm 
dcvcloped herc. First, the algorithm could be developed as a sub-routine to be embedded 
in a larger device simulation code. In the accurate modeling of PN junctions, for 
cxample, especially those utilizing heavily doped semiconductor regions (as in Esaki 


diodes), one necds to make repeated calls to a gencral function that evaluates the FD 
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integrals. If such function calls become too numerous, however, one might envision the 
use of our algorithm as follows. One would use the sub-routine to first evaluate the FD 
integral over a range of values. These values would then be stored as a look-up table, to 


which standard numerical interpolation methods would then be applied. 


Zh 





APPENDIX A. FERMI DIRAC ALGORITHM IN ANSI C 


/*Thesis Program to solve Fermi-Dirac Intergrals with the method of Levine*/ 
/*Lt. Jeremy Thompson*/ 

/*Compiler: Boreland Turbo C++*/ 

/*File Name: fd.c*/ 

#include<stdio.h> 

#include<math.h> 

#include<stdlib.h> 

#define MAX 15 

double inpvar(double); 

void print_results(double, double); 

double eta(double); /* In accordance with (IAW) Eq. 16*/ 
double eq11(double, double); 

double hypo(double,double,double); /* IAW Eq. 22*/ 
double gamma(double); 

double eq18(double,double); 

void fermdir(double); 

double contfrac(double,double); pean 33 7) 
double eq32(double,double); 

double fact(double); 

double levin(); 

double bico(int,int); PoRer 327 
double factin(int n); 

double gammIn(float xx); 


double TIMAX+1],t[MAX+1}; 


void main() 
{double tau; 
tau=inpvar(tau); 
fermdir(tau); 
return: } 


double inpvar(double tau) 
{printf("\n What is the value of tau ?\n"); 
scanf("%olf",&tau); 
return tau; } 


void print_results(double x, double fx) 
{printf("%lf \t%34.321f \n",x,fx); 
retum; } 


void fermdir(double tau) 
{double x,fx,a,b; 
x=-4.0; 
while(x<4.1) 
{if(x<0.0) 
{fx=eq11(x,tau); } 
if(fabs(x)<.01) 
{fx=eta(taut+1.0); } 
if(x>=.1) 
{a=fx=eq 18(tau,x):; 
b=eq32(tau,x); 
fx=fx-b; 
fx=1+(taut+] )*fx; 


Zo 


fx=fx*pow(x,taut+1); 


fx=fx/gamma(tau+2); } 
print_results(x,fx); 
x=x+0.1;} 
return; } 


double eq11(double x, double tau) 
{double y; 
int n,k; 
T[0]=t{0]=0.0; 
n=k=1; 
y=pow(-1,n+1)*exp(n*x)/(pow(n,tau+1)); 
while(k<=MAX) 
{t{k]=pow(-1,n+1)*exp(n*x)/(pow(n,tau+1)); 
T[k]J=T[k-1]+t[k], 
n=n+1; 
k=k+1;} 
y=levin(); 
retum y; $ 


double hypo(double a,double c,double z) 
{long double al,c1,z1,temp,kfact; 


double m; 
int k; 
k=0;al=a;cl=c;kfact=1;temp=0; 
m=1;z1=z; 
while(fabs(m-temp)>1e-14) 
{k=k+1; 
kfact=kfact*k; 
temp=m; 
if(k<=1) 
{m=m-+(a/c*z1/kfact); } 
else 
{a=a*(al+k-1); 
c=c*(cl+k-1); 
m=m+(a/c*z1/kfact); 
} 
z1=z1*z;} 
retum m;} 


double gamma (double tau) 
{double y,temptau,gammatau,temp; 
y=tau; 
temptau=1.0000; 
while(tau>2) 
{tau=tau- 1; 
temptau=temptau* tau, 


s 
while(tau<-2) 
{tau=taut]; 
tcinptau=tcmptau/(tau-1); 


' 
§ 


gamimatau=1 .O0000000 *pow(tau, 1); /* Ret. 2*/ 
gammatau = gammatau + (0.5772 1566490 15329*pow(tau,2)); 
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gammatau = gammatau + (-.65587807 15202538 *pow(tau,3)); 
gammatau = gammatau + (-.0420026350340952*pow(tau,4)); 
gammatau = gammatau + (0.16653861 13822915*pow/(tau,5)):; 
gammatau = gammatau + (-.0421977345555443*pow(tau,6)); 
gammatau = gammatau + (-.0096219715278770*pow/(tau, 7)); 
gammatau = gammatau + (0.0072189432466630*pow(tau,8)); 
gammatau = gammatau + (-.0011651675918591*pow(tau,9)); 
gammatau = gammatau + (-.0002152416741149*pow/(tau, 10)); 
gammatau = gammatau + (0.0001280502823882*pow/(tau, 11)); 
gammatau = gammatau + (-.0000201348547807*pow/(tau, 12)); 
gammatau = gammatau + (-.0000012504934821*pow(tau, 13)); 
gammatau = gammatau + (0.0000011330272320*pow/(tau, 14)); 
gammatau = gammatau + (-.0000002056338417*pow(tau, 15)); 
gammatau = gammatau + (0.0000000061160950*pow/(tau, 16)):; 
gammatau = gammatau + (0.0000000050020075*pow/(tau, 17)); 
gammatau = gammatau + (-.0000000011812746*pow(tau, 18)); 
gammatau = gammatau + (0.000000000 1043427*pow/(tau, 19)); 
gammatau = gammatau + (0.0000000000077823*pow(tau,20)); 
gammatau = gammatau + (-.0000000000036968*pow(tau,21)); 
gammatau = gammatau + (0.0000000000005 100*pow(tau,22)); 
gammatau = gammatau + (-.0000000000000206*pow(tau,23)); 
gammatau = gammatau + (-.0000000000000054 *pow(tau,24)); 
gammatau = gammatau + (0.0000000000000014*pow(tau,25)); 
gammatau = gammatau + (0.0000000000000001*pow(tau,26)); 
tau=y; 

=(1/gammatau)*temptau; 
return y; } 


double eq18 (double tau, double x) 
{double a; 
int k; T[O]=t{0]=0.0; 
k=1; 
while(k<=MAX) 
{t{k]=pow(-1,k)*contfrac(2+tau,k*x); 
T[k]=T[k-1]+t[k];k=k+1;} 
a=levin(); 
return fabs(a); } 


double contfrac(double tau, double z) 
{double a[1000],b[ 1000], value,old; 
int i,k,index; 
index=0; 

tau=tau-2; 
k=0; 
while(k<1000) 

{a{k]=b[k]=0;k++; } 

k=0; 
a{1]=1;b[1 J=z; 
for(i=2;1<=999;++1) 

{if(k<1) 
{afi]=index-tau; 
bfiJ=1; 
index=index+1; 
Kee 


a4 


else 
{a[i]=index; 
b[i]=z; 
k=0;} 


j 
i=1;old=a[i]/b[1]; 
i=2;value=a[i]/b[i]; 
i=i-1;value=a[i]/(b[i]+value); 
ep 
while(fabs(value-old)>1e-16) 
{old=value; 
i=i+1; 
index=1; 
value=a[i]/o[i]; 
while(index>1) 
{index=index-1; 
value=a[index]/(b[index]+value); } 
j 
return value; } 


double eq32(double tau,double x) 
{double value; 
int k; k=1; T[0]=t[0]=0.0; 
while(k<=MAX) 
{t{k]=pow(-1,k+1)*exp(-k*x)*hypo(1+tau,2+tau,k*x); 
T[k]=T [k-1]+t[k];k=k+1;} 
value=levin(); 
value=value/(taut+ 1); 
return value; } 


double levin() 
{double PIMAX+1],Q[MAX+1]; 
int 1;1=1; 
while(l<=MAX) 
{P[IJ=P[1-1]+pow(-1,)*powd,MAX-1)*T[I]/t[l]*bico(MAX.)); 
Qf1J=Q{[I-1]+pow(-1,D*powd,MAX-1)/t[l]*bico(MAX,]1); 
l=]+1;} 
retum (P[MAX]/Q[MAX]); } 


double fact(double k) 
{double kfact; 
kfact=1; 
while(k>1) 
{kfact=kfact*k;k--;} 

return kfact; } 


double cta(double s) 
{double temp; 
int m;m=1; 
while(m<=MAX) 
{T[m]=T[m-1 ]+pow(-1,m+1)/pow(m,s); 
t[{ m}]=pow(-1,m+1)/pow(m,s); 
m=m-+1;} 
temp=levin(); 
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return temp; } 


double bico(int n, int k) jo rxet. 137/ 
{double factln(int n); 
return floor(0.5+exp(factln(n)-factIn(k)-factln(n-k))); 
} 

double factln(int n) jowvet. 137/ 
{double gammlIn(float xx); 


static float a[101]; 

if (n <= 1) return 0.0; 

if (n <= 100) return a[n] ? a{n] : (a[nJ=gammln(n+1.0)); 
else return gammin(n+1.0); 

} 


double gammln(float xx) (Ret l37/ 

{double x,y,tmp,ser; 

static double cof[6]={ 76. 18009172947 146,-86.5053203294 1677, 
24.01409824083091,-1.231739572450155, 
0.1208650973866 179e-2,-0.5395239384953e-5}; 

int j; 

Y=X=XX; 

tmp=x+5.5; 

tmp -= (x+0.5)*log(tmp); 

ser=1.000000000190015; 

for (j=0:j<=5;j++) ser += coffj]/+ty; 

retum -tmp+log(2.50662827463 10005*ser/x); } 
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APPENDIX B. BLAKEMORE'S TABULATED VALUES 
FOR HALF INTEGER VALUES OF f.(x) 
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APPENDIX C. TABULATED COMPARISON TO BLAKEMORE'S DATA 
FOR HALF INTEGER VALUES OF 1 


Tau = -0.5 


Per Cent Error 
0.012% 
0.68317 0,0000% 


[Per Cent Error 
0.003% 

0.68317 [010000% 

3.8 03 
36 0 
35 | 0.02956 [0.029568 | -0.0271% [0.6 | 0.85074 [0.850744 | -0.0005% 
33 | 0.03595 [0.035949] 0,0028% | 0.8 | 0.93826 | 0.938262 | -0,0002% 
3.2 | 0.03962 [0.039625] _-0.0126% | 0.9 | 0.98235 [0.982555 | 0.005% 
3.0 | 0.0481 [0.048103] -0.0062% | 1.1 [ 1.0717 _[1.071666| 0.032% 
29 
27 Id 1.203192 
6 
Is 
2 2.0 
19 [0.13546 [0.135461] -0.0007% | 22 
18 | 0.14826 | 0.14826 | 0.000% | 2.3 | 1.5872 [1.587248| -0.0030% 
2.4 
16 | 0.1772 [0.177121 | —-0.0006% | 2.5 | 1.6663 | 1.666313 | _-0.0008% 
2.6 
1.0 | 0.29402 [0.294028] _-0.0027% [3.1 | 1.8891 _| 1.889147] -0.0025% 
[0.9 | 0.31845 [0.318452] -0.0006% | 3.2 | 1.9242_[ 1.924237] -0.0019% 
[0.8 | 0.34438 5.3 
0.6 0.40077 [osoo7| —_-0.0002% | 3.3 
0.4 | 0.46318 [0.463176] 0.0009% | 3.7 
02 59 
0.0 | 0.6049 [o.cossoa] 0.000 {| [| | | 


Table C1: Comparison to Blakemore's data to the Levinized values in the interval -4 <x <4 for t=-1/2. 
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Tau = 0.5 
Per Cent Error Levine 
[=4.0 [| 0.018199 [0.018198] —0.0055% | 0.1 | 0.82756 [0.827557 
02 
x 
3.5 | 0.02988 | 0.02988 | 0.0000% | 0.6 | 1.2003 | 1.200251 0.0041% 
r-3.3 | 0.036412 | 0.036412 | _0.0000% [08] 1.3791 _| 1.379123] -0.0017% 
[3.2 | 0.040187 [0.040187] _0.0000% | 0.9 | 1.4752 | 1.475161] 0.026% _ 
r=3.1 | 0.044349 [0.044349 | 0.0000% [1.0 | 1.5736 | 1.575641 | -0.0026% _ 
[3.0 | 0.088933 [0.048934] -0.0020% | 1.1 | 1.6806 | 1.680576] 0.014% _ 
13 
13 
7 
2.2 [0.10671 [0.106707] 0.028% | 1.9 | 2.6794 | 2.679391] 0.003% 
a 
=1.9 [0.14225 [0.142247 0.0021% | 2.2 | 3.1249 | 3.124871] _0.0009% 
23 
r=1.6 [ 0.18889 [0.188887] 0.016% [2.3 
“4 [ 0.22759 2.1 
2.9 
03078 
[0.9 [0.35841 [0.358407] 0.008% [3.2 | 4.8653] 4.865358] -0.0012% _ 
(0.6 | 0.46595 [0.465949] 0.002% [3.5] 5.458 | 5.458044 | -0.0008% _ 
3.8 
0.651606 
4.0 
0.0 | 0.76515 [0.765147] 0.008% | [| 


Table C2: Comparison to Blakemore's data to the Levinized values in the interval -4 < x < 4 for t = 1/2. 


oe 
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Tau = 3/2 

[4.0 | 0.018256 | 0.018257 | -0.0035% [O01] 0.9468 | 0.946803 
[=3.5 | 0.030037 | 0.030038 | -0.0033% [0.6 | 14494 | 1.449432 |" -0.0022% _ 
3:3 | 0.036645 | 0.036646 | -0.0027% [0.8 | 17071 [1.707078 | 0.0013% _ 
[-3.2 | 0.040473 | 0.040473 | 0.0000% [0.9 | 1.8497 | 1849755 |" -0.0030% | 
[-3.1 | 0.044696 | 0.044696 | 0.000% [1.0 | 2.0023 | 2.002258 | 0.0021% _ 
0.0010% 

-0.0014% 

0.0017% 

3.1486 | 3.148555 | 0.0014% 


3.3823 0.002% 
[2.3 | 0.098544 | 0.098544 | 0,0000% | 1.8 | 3.6294 0,0001% 
5.890294 | _0.0002% 
19] 0.14581 | 0.145814 | -0.0027% [2.2 | 4.76 | 4.759999 | 0.000% 
1.6 [0.19517 | 0.195172 | -0.0010% [2.5 | 5.7689_| 5.768879 | 0.004% 
0.004% 
0.002% 


0.002% 
0.005% 
PLT] 0.31533_| 0.315327 | 0.010% [3.0 | 7.7886 _| 7.788611 | 0.001% 
1.0 [| 0.34667 | 0.346675 | -0.0014% ] 3.1 | 8.2467 | 8246693 | 0.000% 
[0.9 | 0.38096 _[ 0.380965 | -0.0013% [3.2] 8.7237 _| 8723665 | 0.004% 
0.8 | 0484s | 0.41844 | 0.000% [3.3] 9.2199 | 9.21988 | 0.002% 
9.7351 0.002% 
0.6] 0.504 | 0.504001 | -0.0002% [3.5] 10.271 _|10.271411 | -0.0040% 


0.0037% 
0,0001% 
0.040% 
0.0017% 
[0.1 | 0.79365 _| 0.793647 | 0.000% | 4.0 | 13.26 | 13.260488 | 0.037% 
roof 08672 | 08672 [| 0.000% | [| 


Table C3: Comparison to Blakemore's data to the Levinized values in the interval -4 <x < 4 for t = 3/2. 
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lI 
Lh 
— 
ho 


Be23 32 3682 19 -0.0006% 
Cpe) ey) 3.515476 0.0007% 
3.8192 3.819185 0.0004% 


0.066813 1.4 
2.6 | 0.073795 13 

0.081501 
2.4 | 0.090006 


to 
Ln 


E is 


ee 


4.1456 4.145618 -0.0004% 


Levine | Per Cent Error lakemore| Levine 
1.0182 
3.9 L.U71 | 1.117129 | -0.0026% 
1.225 _| 1.224998 
3.7 1.3425 
1.4704 
“3.5 | 0.030118 | 0.030117 | __0.0033% | 0.6 |_ 1.6095 
1.7606 
| -3.3 | 0.036764 | 0.036764 | 0.0000% | 0.8 | _ 1.9246 
3.2 | 0.040617 | 0.040617 | 0.0000% | 0.9 | 2.1023 
2.2948 
2.5031 
1.2 | 2.7282 
2.9712 
| -2.7 | 0.066813 | 0. . __-0.0006% 

9.073795 | | 0.0007% 

| 0.081501 | | 0.0004% 

[0.090006 [-0.0004% | 
44961 
22 4.8719 
5.2746 
2.0 5.7053 
6.1662 
-1.8 | 0.16297 2.3] 6.658 | 6.658043 
-16 | 0.19846 [0.198458 | 0.0010% [2.5] 7.7419_[ 7.741875 | 0.0003% | 
14 
12 
1.0 
09 
0.635828 
0.3 | 0.69923 18.944 
01 21-469 
0.0 [0.92755 [0.927554[ -o00% [| [ | | 


Table C4: Comparison to Blakemore's data to the Levinized values in the interval -4 < x < 4 for t = 5/2. 
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aw 2 

lakemore er Cent Error 
1.0587] 1.058705] -0.0005% 
11654] 11654] 0.000% 
0.022349 1.2824] 1.282429] -0.0023% 
[4107] 1410721[ _-0.0015% 
15313] 1551278] 0.014% 
(=3.5[ 0.030158] 0.030157[0.0033% | 0.6[ __1.7052| 1.705177|_0.0013% 
18736| 1.873578| 0.0012% 
[-3:3] 0.036824[ 0.036823] 0.027% | 0.8] 2.0577] 2.057724 -0.0012% 
(=3.2[ 0.04069] 0.040689] 0.0025% | 0.9| 2.2589] 2.258948] -0.0021% 
10 


0.044961 0.04496] 0.0022% 1.0 
-3.0| 0.049679) 0.049678] 0.0020% 1.1 


2.4787| 2.478679} 0.0008% 
2.7184] 2.71844} -0.0015% 


0.054891] 0.054891 2.9799 | 2.979861} 0.0013% 
3.2647] 3.264676] 0.007% 
0.067009] 0.067008 3.5747| 3.574733] -0.0009% 
3.912 3.911994] 0.0002% 
[25] 0.08mi] o0si7i| 0.000% | 16] 4.2786] 4.278543| 0.013% 
4.6766] 4.676589] 0.0002% 
5.1085] 5.108468] 0.006% 
5.5767| 5.576653] _ 0.008% 
6.0838] 6.083752] 0.008% 
6.6325 
[-19] 0.1486 [0.148603] -0.0020% [2.2] 7.2258 7.225849 -0.0007% _ 
24 
2.6 
28{ 11.903] 1190345] -0.0038% _ 
TI] 0.32822[ 0.328216] 0.0012% [3.0] 13.976] 13.97612[ -0.0009% 
[-10f0.36222|0.362222[ -0.0006% [3.1] 15.127] 15.12710[-0.0007% _ | 
[09] 0.3997] 0.399697 0.008% | 3.2[ 16.36] 16.36058]  -0.0036% _| 
[-0.8] 0.44098{ 0.440984] -0.0009% [3.3] 17.681] 17.68130[-0.0017% 
| -0.6] 0.53653] 0.53653 20.60455| 0.0021% | 
22176 0.018% | 
0,65229| 0.65229| 0.000% | 3.7[ 23.939] 23.93893] 0.003% 
02 39[ __27.73| 27.72975| _0.0009% 
[01] 0.87294] 0.872939] 0.0001% —[ 40[ 29.812 29.81139[ 0.020% 
oo] 0.96148] 0.961484, 0.000% [ [| | 


Table C5: Comparison to Blakemore's data to the Levinized values in the interval -4< x < 4 for t = 7/2. 
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